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We calculate thermodynamic properties of soft-core lattice bosons with on-site n-body interac- 
tions using up to twelfth and tenth order strong coupling expansion in one and two dimensional 
cubic lattices at zero temperature. Using linked cluster techniques, we show that it is possible to 
exactly renormalize the two-body interactions for quasiparticle excitations and ground-state energy 
by resumming the three and four body terms in the system, which suggests that all higher-body 
on-site interactions may be exactly and perturbatively resummed into the two-body terms. Such a 
procedure is applicable to a broad range of systems analyzable by linked cluster expansions, ranging 
from perturbative quantum chromodynamics to spin models. Universality at various three-body 
interaction strengths for the two dimensional boson Hubbard model is checked numerically. 
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The first calculation of the effect of three-body inter- 
actions in lattice bosons [l| for liquid He 4 and solid He 
revealed its negligible effect on its ground state energy. 
However, it was recently suggested [2j that, firstly, three- 
body interactions in cold polar molecules could be nat- 
urally modelled by Hubbard Hamiltonians with nearest- 
neighbour three-body terms and, secondly, there might 
arise plausible new phases in experimental setups of de- 
generate quantum molecular gases. Shortly thereafter, a 
decoupling mean-field (MF) approximation was used to 
investigate the critical properties of a boson Hubbard 
model with on-site three-body interactions That 
such on-site terms could effectively arise in two-body 
collisions of atoms confined to optical lattices was only 
subsequently justified We will use a strong coupling 
expansion - with no finite size effects in the thermody- 
namic limit - of the on-site three-body interacting boson 
Hubbard model, in addition to checking the universality 
hypothesis at XY critical points for various three-body 
strengths. A procedure for incorporating higher body in- 
teractions by renormalizing the two-body problem will 
also be described herein. 

Consider a system of two-body and higher-body inter- 
acting bosons on the one dimensional chain and two di- 
mensional square lattices described by the Bosc-Hubbard 
Hamiltonian 

H = -tJ2 ib\b J + h.c) + - X ) 
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where the b\ and b i are bosonic creation and annihila- 
tion operators, n,; = is the number operator, the 
hopping-terms t are between nearest neighbors, and the 
system consists of a single species of soft-core bosons. 
The local energy term U contributes to a repulsive on- 



site interaction between bosons, Uk > is the strength 
of fc-body on-site interaction terms and p, is the chemical 
potential. The onsite term U will be the energy scale of 
choice in this letter. 

The Hamiltonian in Eq. (JTJ), when represented in the 
form H = Ho — XHi with Hi being the hopping terms 
of strength A = t/U and Ho being the rest of Eq. (JTJ), is 
amenable to linked cluster expansions 0, Q in the param- 
eter A. Evaluation of physical properties e.g. energy are 
performed via Rayleigh-Schrodinger perturbation theory 
[7] and the linked cluster expansion. Excited states can 
be obtained using a similar procedure through Gelfand's 
similarity transformation 

Consider first the p = 1 Mott insulating lobe in the one 
dimensional chain. For a twelfth order bond-expansion, 
there are 13 distinct topological graphs (clusters) that 
can be embedded on the infinite chain: approximately 
2.5 million states contribute to the full Hilbert space 
with a maximum of 13 states in the lowest degenerate 
manifold. From MF calculations 0, Q , it was predicted 
that the first Mott lobe should not change in structure; 
this was later systematically corrected by densit y m atrix 
renormalization group (DMRG) calculations Us- 
ing Gelfand's similarity transformations to construct the 
particle and hole excitations, we identify the disappear- 
ance of the excitation gap as defining the second-order 
transition contours of the lobe: we have thus evaluated a 
finite series for the gap up to twelfth order. To illustrate, 
for r 3 = ^ = 1, the first four out of twelve calculated 
terms in the Mott gap, Si(X,k = 0), is given by 
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where k is the lattice momentum. In Fig. [T] we show 
particle and hole contours obtained by multiple precision 
Pade approximation of twelfth order series and compared 
to a previously published DMRG solution for = 7. The 
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FIG. 1: First Mott lobe in the ground state phase diagram 
of the one dimensional boson Hubbard model with various 
r-i = ^ ratios; the last Pade approximant [6/6] or [5/5] of 
the twelfth order gap series for r$ — 7 is compared with a 
DMRG solution [m/n] denotes an m th order numerator 
and an n order denominator. The Kosterlitz-Thouless point 
shifts upwards and rightwards as 7-3 is increased. 



location of the critical point (Kosterlitz-Thouless tran- 
sition) shifts upwards and rightwards in the phase dia- 
gram indicating an increase in the size of the first lobe 
and a partial restoration of particle-hole symmetry as 
the semi-hardcore condition (ra = 00) is reached. We 
have verified this tendency with = 0, 1, 7, 100. For the 
hardcore condition, the critical \ij~U (U now being the 
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nearest neighbour repulsion) equals exactly 1 
the particle-hole symmetry completely restored. 
We note that in an n-th order bond expansion for the lin- 
ear chain, we need calculate the n-th order particle and 
hole contributions only up to graphs with n— 1 bonds: the 
effective Hamiltonian for the last cluster may be calcu- 
lated with very little effort because, within the degener- 
ate subspace of this graph, the only contributing process 
will be the one which transfers the excitation from one 
end of the chain to the other. These matrix elements will 
simply be given by the negative of the Schroeder num- 
bers S = {1,2,6, 22, 90 • • • } with the generating function 
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That is, for a cluster with n-bonds the n order ef- 
fective Hamiltonian has its non-zero elements given by 
H(0,n) = —S n -i, for I > 2. This is independent of r% 
and the type of excitation. 

In two dimensions there are 680 topologically distinct 
clusters contributing at tenth order for a bond-expansion. 
Here the Mott gap closes as 8 cx (t c — t) zl ' [13| for 
t — t c -C 1, assuming the universality of the XY model: 
t c is the value of the hopping element at the multicriti- 
cal point where particle-hole symmetry is restored (here 



TABLE I: Critical points and exponents for the two dimen- 
sional square lattice boson Hubbard model at various three- 
body interactions 7-3 using roots extrapolation (E) and Pade 
approximations (P). See text and Ref. |(| for the exact pro- 
cedure. At t — tc, z — 1 Hal. 
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597.4 ± 0.4 
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690 


1 


603.8 ± 0.8 


604.69 ± 0.06 


692.3 ± 0.6 


10 


616.7 ± 0.8 


617.39 ± 0.05 


695.4 ± 0.4 


100 


621.4 ± 0.8 


621.98 ± 0.06 


696.5 ± 0.5 
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z = 1), z and v are the dynamical and coherence length 
critical exponents respectively. To investigate the effect 
of three-body terms in 2D, we have calculated tenth or- 
der series for the p = 1 Mott gaps for r$ = 1,10,100. 
From these series', t c and zv may be extracted by pro- 
ceeding, mutatis mutandis, as outlined in previous scaling 
analysis 0] : (a) by linearly extrapolating the roots of the 
truncated gap series' from, say, fourth to tenth order and 
(b) by Pade approximating the gap series' to mimic the 
expected 6 behaviour mentioned above. The results of 
the higher approximants ([4/4], [4/5], [4/6], [5/4], [5/5]) 
and linear extrapolation arc tabulated in Table |T] for four 
r3 values: it must be noted that large r% values may be 
attained, as suggested by Johnson et. al, using Feshbach 
resonances and tuning the lattice potential. As can be 
noted from the table that the change in t c upon increas- 
ing r$, and hence the structure of the first lobe, is not as 
substantial for the square lattice as was seen for the one 
dimensional case. From the v values for the four three- 
body interaction strengths, we see that universality does 
indeed seem to hold at the XY point. The corresponding 
classical critical coefficient for the three dimensional XY 
model is v = 0.67155 ± 0.00027 [l4|. 



In general, to incorporate a second variable like U3 
into a Hamiltonian within linked cluster expansions re- 
quires a double-expansion: the first in t/U, the second 
in U3/U. For example, the double-expansion of a quan- 
tity P in perturbing variables A and r to order M and N 
respectively may be symbolically written as 
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Now M and N are finite integers but can one do bet- 
ter? The prescription we adopt is to resum the second 
series and evaluate lim.Ar_ ) . 00 c^" 1 for every i, keeping M 
finite, and is implemented as follows: we first calculate 
the series coefficients for a given observable (like in Eq. (2) 
for a finite number of T3 values. And because the coef- 
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TABLE II: Resummed series coefficients for the particle and hole contours in the one dimensional chain and two dimensional 
square lattice. Coefficients of lower order that are not listed are independent of r 3 . 
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ficients are always rational numbers - by virtue of the 
perturbation theory - it only remains to find a rational 
function approximation to the obtained coefficients. The 
latter step may be easily implemented with Thiele's al- 
gorithm for continued fraction representation (l5| . 
To summarize, the functional dependence of a coefficient 
at a given order on r 3 is to be captured by a rational ap- 
proximation. For example, a single-expansion coefficient 
Cu, for a given i, for some 24 values of r 3 from — > 100 
were evaluated. Thiele's algorithm to find an approxi- 
mating rational function /j(r 3 ) = Cu(r 3 ) would generally 
require as many steps as there are points (here 24) to ter- 
minate and find the best fit; however, we find that in each 
of the evaluated coefficients, the algorithm stops exactly 
after a few steps because the continued fraction expan- 
sions stop. This ensures the exactness of the obtained 
fiirz). With this, the Ci's in Eq.[3]get fully renormalized 
by the resummed cVs. 

The above procedure has been applied to renormalizing 
the scries coefficients of the particle and hole contours in 
the two-body interacting one dimensional chain and two 
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FIG. 2: Fourth order coefficient for the hole contour in the one 
dimensional chain as a function of the three-body interacting 
strength in a log plot. The coefficient at r 3 = passes through 
the function as well. The [1/1] function for is 6 °~^ 3 . 



dimensional square lattice with respect to the three body 
terms. Let the particle and hole contours, for any r 3 , be 
represented as 



±/i c ± (r 3 ) = ^c±(r 3 )A i . 



(4) 



i=0 



The signs (±) refer to the particle and the hole contours 
respectively. For illustrating our method, we sketch in 
Fig.[2]the [1/1] rational function approximation to (r 3 ) 
in the one dimensional chain obtained from the 24 dif- 
ferent values of r 3 . The same analyses were performed 
for the particle coefficients as well and similar conclu- 
sions hold; the resummed coefficients are listed in Ta- 
ble |TT] up to fourth order. For example, in the ID case, 
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Eq. [21 It is worth noting that even with coefficients for 
particle-hole series only up to third order, very reason- 
able estimates (within 10% compared to more accurate 
results 0) for critical properties can be obtained flij ]. 
A similar rcsummation may be readily obtained, exactly 
and without using the above procedure, for the ground 
state energy per site (tp\'H\'(p) - where \ip) is the ground 
state wavefunction constructed order by order for the first 
lobe - in the one dimensional chain for an arbitrary r up 
to fourth order to give 



E}, D = -4A 2 



1 2 1±^A 4 . 
3 + r 3 
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We see from Table [TT] and Eq. [5] that for certain values 
of attractive interactions i.e. r 3 < there is a per- 
turbativc instability of the p — 1 Mott phase coming 
from the divergence of the denominators. This might sig- 
nal the disappearance of the first lobe altogether or the 
appearance of a higher-density and energetically more 
favourable lobe in that region of phase space: quite 
naturally, for attractive bosons, higher density Mott 
phases should stabilize the system and one should ex- 
pand thermodynamic variables perturbatively about this 
more favourable phase. Similar conclusions were in fact 
reached by recent MF and quantum Monte Carlo calcu- 
lations 17[. In the present work, however, the value of 
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the attractive three-body strength that leads to an insta- 
bility at a given perturbative order can be readily read 
off from the resummed coefficients. 

Using the above procedure for the Hamiltonian with four 
body interactions, with r± = Ht-, we find that the two- 



body interactions for the ground state energy densities of 
the linear chain and the square lattice may also be per- 
turbatively renormalized by the four-body terms as given 
by the following: 



1D 2 x4 272 6 4(85r 4 - 7602) 8 2(252109r| + 2870730^ + 6509628r 4 - 9540936) . 10 
E = -4A + 4A + — A + — — — — — — A — — A +••• 



9 81(r 4 + 6) 729(r 4 + 6) 3 

I- 27 4 39 514r| + 6333r| + 28167ri + 531 
r 4 + 3~ ~ " (r 4 + 3) 3 (2r 4 + 3) 



E 2B = _ 8X 2 _ 21 7r 4 + 27 A4 _ 32 514rj + 6333r 3 + 28167rl + 53160r 4 + 35217 , 6 + _ _ _ _ (g) 



Therefore it seems very likely that two-body terms in the 
Bosc-Hubbard model, irrespective of dimension, may be 
perturbatively renormalized by all higher-body on-site 
terms for its thermodynamic and excited properties. An 
interesting question is if such resummability might also 
exist for dynamical properties and for bosonic models 
with intersite interactions. 

In conclusion, we have presented a way of resumming the 
effect of a second perturbing variable to infinite order 
thereby effectively renormalizing the series coefficients of 
the single variable expansions. The procedure is quite 
general and may be applied to renormalizc the second in- 
teraction term in the Hamiltonian in the series expansion 
representation of any thermodynamic quantity. In the 
scenario considered, we have found perturbative renor- 
malizations of the two-body interactions due to three and 
four body terms in calculations of ground state energies 
and quasiparticle excitations, for the one and two dimen- 
sional Bose-Hubbard model. The applicability of the pro- 
cedure is, of course, not restricted to lattice bosons but 
can be extended to other systems that are treated using 
series expansion techniques, ranging from spin models to 
perturbative quantum chromodynamics where the ana- 
lytic continuation of str ong coupling expansions is still 
fraught with problems [181 ]. Additionally, the univer- 
sality hypothesis of the two dimensional Bosc-Hubbard 
model has been checked for various three body interac- 
tion strengths. 

One of us (VKV) thanks the Bonn-Cologne Graduate 
School for support within the Deutsche Forschungge- 
meinschaft's Research Funding. 



[2] H.P. Buechler, A. Micheli, and P. Zoller. Nat. Phys., 
3:726, 2007. 

[3] B 1. Chan, Xiao bin Huang, Su peng Kou, and Yunbo 

Zhang. Phys. Rev. A, 78:043603, 2008. 
[4] P.-R. Johnson, E. Tiesinga, J. V. Porto, and C. J. 

Williams. New Journal of Physics, 11:093022, 2009. 
[5] M.P Gelfand and R.R.P. Singh. Adv. in Phys., 49:93, 

2000. 

[6] N. Elstner and H. Monien. Phys. Rev. B, 59:12184, 1999. 
[7] Gordon Baym. Lectures on Quantum Mechanics. 

Lecture Notes and Supplements in Physics. Ben- 

jamin/Cummings, England, 1969. 
[8] Kezhao Zhou, Zhaoxin Liang, and Zhidong Zhang. Phys. 

Rev. A, 82:013634, 2010. 
[9] J. Silva- Valencia and A.M.C Souza. Phys. Rev. A, 

84:065601, 2011. 
[10] Manpreet Singh, Arya Dhar, Tapan Mishra, R. V. Pai, 

and B. P. Das. Phys. Rev. A, 85:051604(R), 2012. 
[11] C.N. Yang and CP. Yang. Phys. Rev., 151:258, 1966. 
[12] The on-line encyclopedia of integer sequences. 

|http://oeis . org/A006318| 2012. 
[13] Mathew P.A. Fisher, Peter B. Weichman, G. Grinstein, 

and Daniel S. Fisher. Phys. Rev. B, 40:564, 1989. 
[14] Massimo Campostrini, Martin Hasenbusch, Andrea 

Pelissetto, Paolo Rossi, and Ettore Vicari. Phys. Rev. 

B, 63:214503, 2001. 
[15] J. Stoer and R. Bulirsch. Introduction to Numerical Anal- 
ysis. Texts in Applied Mathematics. Springer, 3 edition, 

2002. 

[16] J.K. Freericks and H. Monien. Phys. Rev. B, 53:2691, 
1996. 

[17] A. Safavi-Naini, J. von Stecher, B. Capogrosso-Sansone, 
and Seth T. Rittenhouse. Phys. Rev. Lett., 109:135302, 
2012. 

[18] Jaan Oitmaa, Chris Hamer, and Weihong Zheng. Series 
Expansion Methods for Strongly Interacting Lattice Mod- 
els. Cambridge University Press, Cambridge, 1 edition, 
2006. 



[1] RD. Murphy and J.A. Baker. Phys. Res. A, 3:1037, 1971. 



